#include "Modeling\CUDADataStruct.h"
#include ".\CUDA\cutil\inc\cutil.h"
#include <stdio.h>

#define PI 3.1415926536f


#define USE_SHARED 1

__constant__ My_CUDA_Global* d_globalData;
__constant__ My_CUDA_Sil* d_silData;
__constant__ int* d_silX;
__constant__ float* d_ZStep;
__constant__ float* d_silZ;
//__constant__ float* d_silUVW;
__device__ float* d_silUVW = 0;

typedef unsigned int uint;
uint degX, degY, degZ;


int iDivUp(int a, int b)
{
    return (a % b != 0) ? (a / b + 1) : (a / b);
}

__global__ void cuda_kernel_texture_2d(unsigned char* surface, size_t pitch, int width, int height)
{
    int x = blockIdx.x*blockDim.x + threadIdx.x;
    int y = blockIdx.y*blockDim.y + threadIdx.y;
    //float* pixel;
    unsigned char* pixel;
    float t=0.0f;
       
    // in the case where, due to quantization into grids, we have
    // more threads than pixels, skip the threads which don't 
    // correspond to valid pixels
	if (x >= width || y >= height) return;
	
    // get a pointer to the pixel at (x,y)
    //pixel = (float*)(surface + y*pitch) + 4*x;
    pixel = surface + y*pitch + 4*x;

	// populate it
	float value_x = 0.5f + 0.5f*cos(t + 10.0f*( (2.0f*x)/width  - 1.0f ) );
	float value_y = 0.5f + 0.5f*cos(t + 10.0f*( (2.0f*y)/height - 1.0f ) );
	//pixel[0] = 0.5*pixel[0] + 0.5*pow(value_x, 3.0f); // red
	//pixel[1] = 0.5*pixel[1] + 0.5*pow(value_y, 3.0f); // green
	//pixel[2] = 0.5f + 0.5f*cos(t); // blue
	//pixel[3] = 1; // alpha

	pixel[0] = (0.5f + 0.5f*cos(t))*255; // blue
	pixel[1] = (0.5*pixel[1] + 0.5*pow(value_y, 3.0f))*255; // green
	pixel[2] = (0.5*pixel[0] + 0.5*pow(value_x, 3.0f))*255; // red
	pixel[3] = 255; // alpha
}

extern "C" 
void cuda_texture_2d(void* surface, size_t pitch, int width, int height)
{
    dim3 Dg = dim3(32,32);
    dim3 Db = dim3( (width+Dg.x-1)/Dg.x, (height+Dg.y-1)/Dg.y );
    cuda_kernel_texture_2d<<<Dg,Db>>>( (unsigned char*)surface, pitch, width, height);
}

__device__ void normalize(float* a)
{
	float sum;
	sum=sqrt(a[0]*a[0]+a[1]*a[1]+a[2]*a[2]);
	if(sum > 1.0e-10f)
    {
		a[0]/=sum; a[1]/=sum; a[2]/=sum;
	}
}

__device__ void addambient(MATERIAL *mat, float rgb[])
{
	rgb[0]+=mat->ka*mat->d[0]*mat->a[0];
	rgb[1]+=mat->ka*mat->d[1]*mat->a[1];
	rgb[2]+=mat->ka*mat->d[2]*mat->a[2];
}

__device__ void addPhong(MATERIAL *mat, LIGHT* light, float norm[], float ray[], float rgb[])
{
	float H[3] ={(-light->dir[0]-ray[0])/2.0f, (-light->dir[1]-ray[1])/2.0f, (-light->dir[2]-ray[2])/2.0f};
	normalize(H);
	float phong = norm[0]*H[0]+norm[1]*H[1]+norm[2]*H[2];
	phong *= phong;
    phong *= phong;
    phong *= phong;
    phong *= mat->ks;
	rgb[0]+=phong*light->rgb[0];
	rgb[1]+=phong*light->rgb[1];
	rgb[2]+=phong*light->rgb[2];
}

__device__ void addDiffuse(MATERIAL *mat, LIGHT* light, float diff, float norm[], float rgb[])
{
	rgb[0]+=mat->kd*diff*light->rgb[0];
	rgb[1]+=mat->kd*diff*light->rgb[1];
	rgb[2]+=mat->kd*diff*light->rgb[2];
}

__device__ void GetPlanesUVW(float4* vertplane, float4* horPlane, int x, int y, unsigned int width, unsigned int height, float vpUVW[][4])
{
    float l;
    l=2.0f*x/width-1.0f;
	vertplane->x = (l*vpUVW[0][3]-vpUVW[0][0]);
	vertplane->y = (l*vpUVW[1][3]-vpUVW[1][0]);
	vertplane->z = (l*vpUVW[2][3]-vpUVW[2][0]);
	vertplane->w = (l*vpUVW[3][3]-vpUVW[3][0]);
	
    l=1.0f-2.0f*y/height;
	horPlane->x = (l*vpUVW[0][3]-vpUVW[0][1]);
	horPlane->y = (l*vpUVW[1][3]-vpUVW[1][1]);
	horPlane->z = (l*vpUVW[2][3]-vpUVW[2][1]);
	horPlane->w = (l*vpUVW[3][3]-vpUVW[3][1]);
}

__device__ void GaussEliminationGPU(float A[][3], float b[])
{
	int i,j,k,ik;
	float mik, temp;
	for(k=0; k<3; k++)
	{
		mik=-1.0f;
		for(i=k; i<3; i++)
		{
			if(fabs(A[i][k])>mik)
			{
				mik=fabs(A[i][k]);
				ik=i;
			}
		}
		if (mik<1.0e-10f)
		{
			b[0]=b[1]=b[2]=0.0f;
			return;
		}
		for(j=k; j<3; j++)
		{
			//swap(A[ik][j], A[k][j]);
			temp = A[ik][j];
			A[ik][j]=A[k][j];
			A[k][j]=temp;
		}
		//swap(b[ik], b[k]);
		temp = b[ik];
		b[ik] = b[k];
		b[k] = temp;
		b[k] /= A[k][k];
		for(i=3-1; i>=k; i--)
			A[k][i]/=A[k][k];
		for(i=k+1; i<3; i++)
		{
			b[i] -= A[i][k]*b[k];
			for (j=3-1; j>k; j--)
				A[i][j] -= A[i][k]*A[k][j];
		}
	}
	for(i=3-1; i>=0; i--)
	{
		for (j=i+1; j<3; j++)
			b[i] -= A[i][j]*b[j];
	}
}

__device__ void EvalBezier(float u, float v, float w, int nx, int ny, int nz, float func[][4][4], float f[][4][4], float* re)
{
    int i,j,k,l;
	float cu=1.0f-u, cv=1.0f-v, cw=1.0f-w;
    
    for(k=0;k<=nz;k++)
    for(j=0;j<=ny;j++)
    for(i=0;i<=nx;i++)
        f[i][j][k]=func[i][j][k];

	for(j=0;j<=ny;j++)
		for(k=0;k<=nz;k++) 
			for(l=1;l<=nx;l++)
				for(i=0;i<=nx-l;i++)
					f[i][j][k]=cu*f[i][j][k]+u*f[i+1][j][k];

	for(k=0;k<=nz;k++) 
		for(l=1;l<=ny;l++)
			for(j=0;j<=ny-l;j++)
				f[0][j][k]=cv*f[0][j][k]+v*f[0][j+1][k];


	for(l=1;l<=nz;l++)
		for(k=0;k<=nz-l;k++)
			f[0][0][k]=cw*f[0][0][k]+w*f[0][0][k+1];

	*re=f[0][0][0];    
}

__device__ void GetBernstein(float Bern[], int deg, float u)
{
	float cu = 1.0f-u;
	switch (deg)
	{
	case 2:
		Bern[0]=cu*cu; Bern[1]=2*u*cu; Bern[2]=u*u; break;
	case 1:
		Bern[0]=cu; Bern[1]=u; break;
	case 3:
		Bern[0]=cu*cu*cu; Bern[1]=3*u*cu*cu; Bern[2]=3*u*u*cu; Bern[3]=u*u*u; break;
	}
}

__device__ void EvalBezierBlossom(float u, float v, float w, int nx, int ny, int nz, float func[][4][4], float* re)
{
	float nn[4];
	GetBernstein(nn, nz, w);

	float f1[4][4];
	float f2[4];
	for (int i=0; i<=nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += func[i][j][k]*nn[k];
	}

	GetBernstein(nn, ny, v);

	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nn[j];
	}

	GetBernstein(nn, nx, u);

	(*re) = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		(*re) += f2[i]*nn[i];
	}
}

__device__ void GetPowerBase(float Pow[], int deg, float u)
{
    Pow[0] = 1.0f;
    Pow[1] = u;
    Pow[2] = u*u;
    Pow[3] = Pow[2]*u;
}


__device__ void EvalBezierPower(float u, float v, float w, int nx, int ny, int nz, float func[][4][4], float* re)
{
	float nu[4], nv[4], nw[4];
    GetPowerBase(nu, nx, u);
	GetPowerBase(nv, ny, v);
	GetPowerBase(nw, nz, w);

	float f1[4][4];
	float f2[4];
	for (int i=0; i<=nx; i++)
	for (int j=0; j<=ny; j++)
	{
		f1[i][j] = 0.0f;
		for (int k=0; k<=nz; k++)
			f1[i][j] += func[i][j][k]*nw[k];
	}
	for (int i=0; i<=nx; i++)
	{
		f2[i] = 0.0f;
		for (int j=0; j<=ny; j++)
			f2[i] += f1[i][j]*nv[j];
	}

	(*re) = 0.0f;
	for (int i=0; i<= nx; i++)
	{
		(*re) += f2[i]*nu[i];
	}
}


__device__ void usrfunGPU(float point[], float fValue[], float jacobi[][3], float func[][4][4], float4 scanplane, float4 vertplane, float jacobiF[][4][4][4], int dx, int dy, int dz)
{
    float u=point[0], v=point[1], w=point[2];
    
	EvalBezierBlossom(u, v, w, dx, dy, dz, func, &fValue[0]);
	fValue[1] = scanplane.x*u + scanplane.y*v + scanplane.z*w + scanplane.w;
	fValue[2] = vertplane.x*u + vertplane.y*v + vertplane.z*w + vertplane.w;
	
	//f=jacobiF[0];
	EvalBezierBlossom(u, v, w, dx-1, dy, dz, jacobiF[0], &jacobi[0][0]);
	//f=jacobiF[1];
	EvalBezierBlossom(u, v, w, dx, dy-1, dz, jacobiF[1], &jacobi[0][1]);
	//f=jacobiF[2];
	EvalBezierBlossom(u, v, w, dx, dy, dz-1, jacobiF[2], &jacobi[0][2]);

	jacobi[1][0] = scanplane.x;
	jacobi[1][1] = scanplane.y;
	jacobi[1][2] = scanplane.z;

	jacobi[2][0] = vertplane.x;
	jacobi[2][1] = vertplane.y;
	jacobi[2][2] = vertplane.z;
}

__device__ void SilTraceInterval(int x, float3* point, float ZStepXX, bool* OutDomain, float4 scanplane, float4 vertplane, float func[][4][4], float jacobiF[][4][4][4], float vpUVW[][4], My_CUDA_Global* globalData)
{
	float pointNew[3] = {point->x, point->y, point->z};
	float fValue[3]={0.0, 0.0, 0.0};
	float jacobi[3][3] = {{0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}};
	float errF, deltaPoint;
    
    float InitDirect[3];
	float TriPlane[3][3];
	float T1P = vpUVW[0][0]*pointNew[0] + vpUVW[1][0]*pointNew[1] + vpUVW[2][0]*pointNew[2] + vpUVW[3][0];
	float T3P = vpUVW[0][2]*pointNew[0] + vpUVW[1][2]*pointNew[1] + vpUVW[2][2]*pointNew[2] + vpUVW[3][2];
	float T4P = vpUVW[0][3]*pointNew[0] + vpUVW[1][3]*pointNew[1] + vpUVW[2][3]*pointNew[2] + vpUVW[3][3];

	TriPlane[0][0] = globalData->renderWidth*(vpUVW[0][0]+vpUVW[0][3]) - 2*x*vpUVW[0][3];
	TriPlane[0][1] = globalData->renderWidth*(vpUVW[1][0]+vpUVW[1][3]) - 2*x*vpUVW[1][3];
	TriPlane[0][2] = globalData->renderWidth*(vpUVW[2][0]+vpUVW[2][3]) - 2*x*vpUVW[2][3];
	InitDirect[0] = 2*x*T4P - globalData->renderWidth*(T1P+T4P);

	TriPlane[1][0] = vpUVW[0][2]-vpUVW[0][3]*ZStepXX;
	TriPlane[1][1] = vpUVW[1][2]-vpUVW[1][3]*ZStepXX;
	TriPlane[1][2] = vpUVW[2][2]-vpUVW[2][3]*ZStepXX;
	InitDirect[1] = T4P*ZStepXX-T3P;

	TriPlane[2][0] = scanplane.x;
	TriPlane[2][1] = scanplane.y;
	TriPlane[2][2] = scanplane.z;
	InitDirect[2] = 0.0f;

	GaussEliminationGPU(TriPlane, InitDirect);
    
	pointNew[0] += InitDirect[0];
	pointNew[1] += InitDirect[1];
	pointNew[2] += InitDirect[2];

   
    for(int its=0; its<4; its++)
	{
		pointNew[0] -= fValue[0];
		pointNew[1] -= fValue[1];
		pointNew[2] -= fValue[2];
		
		//User function supplies function values and Jacobian matrix at point
		usrfunGPU(pointNew, fValue, jacobi, func, scanplane, vertplane, jacobiF, globalData->degX, globalData->degY, globalData->degZ);
		errF = fValue[0]*fValue[0]+fValue[1]*fValue[1]+fValue[2]*fValue[2];
		if (errF <= globalData->errorTol)		break;

		GaussEliminationGPU(jacobi, fValue);
		deltaPoint = fabs(fValue[0])+fabs(fValue[1])+fabs(fValue[2]);
		if (deltaPoint <= globalData->errorXTol)		break;	
	}

    point->x = pointNew[0];
    point->y = pointNew[1];
    point->z = pointNew[2];

    //Compute OutDomain
}

__device__ void PointShading(int x, int y, unsigned char* surface, size_t pitch, float3 point, bool OutDomain, float jacobiF[][4][4][4], float boxX[], float boxY[], float boxZ[], My_CUDA_Global* globalData)
{
    if (!OutDomain)
	{
        // get a pointer to the pixel at (x,y)
        float* pixel;
        pixel = (float*)(surface + y*pitch) + 4*x;
        float rgb[3]={0.0f, 0.0f, 0.0f};
        float u=point.x, v=point.y, w=point.z;
		float norm[3] = {0.0f, 0.0f, 0.0f};
		EvalBezierBlossom(u, v, w, globalData->degX-1, globalData->degY, globalData->degZ, jacobiF[0], &norm[0]);
		EvalBezierBlossom(u, v, w, globalData->degX, globalData->degY-1, globalData->degZ, jacobiF[1], &norm[1]);
		EvalBezierBlossom(u, v, w, globalData->degX, globalData->degY, globalData->degZ-1, jacobiF[2], &norm[2]);
        
		norm[0] /= (boxX[1]-boxX[0]);
		norm[1] /= (boxY[1]-boxY[0]);
		norm[2] /= (boxZ[1]-boxZ[0]);
		normalize(norm);
        
        point.x = point.x*(boxX[1]-boxX[0])+boxX[0];
        point.y = point.y*(boxY[1]-boxY[0])+boxY[0];
        point.z = point.z*(boxZ[1]-boxZ[0])+boxZ[0];

		float ray[3]={point.x-globalData->eye[0], point.y-globalData->eye[1], point.z-globalData->eye[2]};
		normalize(ray);

        LIGHT light;
        float diff;
		for(int i=0; i<globalData->nLight; i++)  
		{
            light.dir[0]=globalData->light[i].dir[0];
            light.dir[1]=globalData->light[i].dir[1];
            light.dir[2]=globalData->light[i].dir[2];
            light.rgb[0]=globalData->light[i].rgb[0];
            light.rgb[1]=globalData->light[i].rgb[1];
            light.rgb[2]=globalData->light[i].rgb[2];
            //normalize(light.dir);
	        diff = -(norm[0]*light.dir[0]+norm[1]*light.dir[1]+norm[2]*light.dir[2]);
	        if (diff>0)
            {
                addPhong(&globalData->material, &light, norm, ray, rgb);
                addDiffuse(&globalData->material, &light, diff, norm, rgb);
            }
		}
		addambient(&globalData->material, rgb);
        pixel[0] = rgb[0];
        pixel[1] = rgb[1];
        pixel[2] = rgb[2];
	}
}

__global__ void silRenderKernel(unsigned char* surface, size_t pitch, int degX, int degY, int degZ, My_CUDA_Global* d_globalData, My_CUDA_Sil* d_silData, int* d_silX, float* d_ZStep, float* d_silZ, float* d_silUVW)
{
	unsigned int x, y, z, i;
    
#if USE_SHARED
    // use shared memory to avoid using local
	__shared__ float sharedFunc[4][4][4];
	__shared__ float sharedJacobiF[3][4][4][4];
	
	//x = threadIdx.x/4 + threadIdx.y/4*2;
	//y = threadIdx.y%4;
	//z = threadIdx.x%4;

	i = threadIdx.y*blockDim.x + threadIdx.x;
	x = i & 0x00000003;
	y = (i>>2) & 0x00000003;
	z = (i>>4) & 0x00000003;
	sharedFunc[x][y][z] = d_silData->func[x][y][z];
	
     __syncthreads();
     
    if (x<degX)
        sharedJacobiF[0][x][y][z] = degX*(sharedFunc[x+1][y][z]-sharedFunc[x][y][z]);
    else
		sharedJacobiF[0][x][y][z] = 0;
    if (y<degY)
		sharedJacobiF[1][x][y][z] = degY*(sharedFunc[x][y+1][z]-sharedFunc[x][y][z]);
    else
		sharedJacobiF[1][x][y][z] = 0;		
	if (z<degZ)
		sharedJacobiF[2][x][y][z] = degZ*(sharedFunc[x][y][z+1]-sharedFunc[x][y][z]);
    else
		sharedJacobiF[2][x][y][z] = 0;
		
	//Even more slow!	
    //sharedJacobiF[0][x][y][z] = degX*(sharedFunc[max(x+1,degX)][y][z]-sharedFunc[x][y][z]);
	//sharedJacobiF[1][x][y][z] = degY*(sharedFunc[x][max(y+1,degY)][z]-sharedFunc[x][y][z]);
	//sharedJacobiF[2][x][y][z] = degZ*(sharedFunc[x][y][max(z+1,degZ)]-sharedFunc[x][y][z]);

    
/*   
    sharedJacobiF[0][x][y][z] = d_silData->jacobiF[0][x][y][z];
    sharedJacobiF[1][x][y][z] = d_silData->jacobiF[1][x][y][z];
    sharedJacobiF[2][x][y][z] = d_silData->jacobiF[2][x][y][z];
*/    
    __syncthreads();
#endif

    x = blockIdx.x*blockDim.x + threadIdx.x + d_silData->minV;
    y = blockIdx.y*blockDim.y + threadIdx.y + d_silData->minSil;
            
    // in the case where, due to quantization into grids, we have
    // more threads than pixels, skip the threads which don't 
    // correspond to valid pixels
	// if (x >= d_globalData->renderWidth || y >= d_globalData->renderHeight) return;    
  
    //skip the pixel out the content of silhouette;
    if (y < d_silData->minSil || y > d_silData->maxSil) return;
    int silN = d_silData->silN;
    if (x <= d_silX[silN*(y-d_silData->minL)] || x >= d_silX[silN*(y-d_silData->minL)+silN-1]) return;

    float3 point;
    float4 verticalPlane, horizonPlane;
    bool OutDomain;
    float ZStepXX, t0;
    for(int interval=silN*(y-d_silData->minL); interval<silN*(y-d_silData->minL)+silN; interval+=2)
    {
        if(x<=d_silX[interval] || x>=d_silX[interval+1]) continue;

        OutDomain=false;
        point.x = d_silUVW[interval*3];
        point.y = d_silUVW[interval*3+1];
        point.z = d_silUVW[interval*3+2];

	    t0 = (0.0f+x-d_silX[interval])/(d_silX[interval+1]-d_silX[interval]);
	    if (t0<0.5f)
	    {
		    t0 *= 2.0f;
		    ZStepXX = (1-t0)*d_silZ[interval] + t0*d_ZStep[interval/2];
	    }
	    else
	    {
		    t0 = 2*(1.0f-t0);
		    ZStepXX = (1-t0)*d_silZ[interval+1] + t0*d_ZStep[interval/2];
	    }
#if USE_SHARED
        GetPlanesUVW(&verticalPlane, &horizonPlane, x, y, d_globalData->renderWidth, d_globalData->renderHeight, d_silData->viewProjUVW);
        SilTraceInterval(x, &point, ZStepXX, &OutDomain, horizonPlane, verticalPlane, sharedFunc, sharedJacobiF, d_silData->viewProjUVW, d_globalData); 
        PointShading(x, y, surface, pitch, point, OutDomain, sharedJacobiF, d_silData->boxX, d_silData->boxY, d_silData->boxZ, d_globalData);
#else
        GetPlanesUVW(&verticalPlane, &horizonPlane, x, y, d_globalData->renderWidth, d_globalData->renderHeight, d_silData->viewProjUVW);
        SilTraceInterval(x, &point, ZStepXX, &OutDomain, horizonPlane, verticalPlane, d_silData->func, d_silData->jacobiF, d_silData->viewProjUVW, d_globalData); 
        PointShading(x, y, surface, pitch, point, OutDomain, d_silData->jacobiF, d_silData->boxX, d_silData->boxY, d_silData->boxZ, d_globalData);
#endif
        break;
    }
     
}

extern "C"
void SilRenderCUDA(void* surface, size_t pitch, int minSil, int maxSil, int minV, int maxV)
{
    dim3 blockNN;
    dim3 threadNN;
    int drawScope = (maxV-minV)*(maxSil-minSil);    

    if(drawScope<=16*64)
    {
        blockNN.x = 4; blockNN.y = 4;
        threadNN.x = iDivUp((maxV-minV), blockNN.x);
        threadNN.y = iDivUp((maxSil-minSil), blockNN.y);        
    }
    else
    {
        threadNN.x = 8; threadNN.y = 8;
        blockNN.x = iDivUp((maxV-minV), threadNN.x);
        blockNN.y = iDivUp((maxSil-minSil), threadNN.y);
    }
	
    
    silRenderKernel<<<blockNN, threadNN>>>((unsigned char*)surface, pitch, degX, degY, degZ, d_globalData, d_silData, d_silX, d_ZStep, d_silZ, d_silUVW);
}

extern "C"
void InitSilRender(My_CUDA_Sil* silData, My_CUDA_Global* globalData)
{
	degX = globalData->degX;
	degY = globalData->degY;
	degZ = globalData->degZ;
    int dSize=silData->silN*(silData->maxL-silData->minL);
	
	CUDA_SAFE_CALL(cudaMalloc((void**)&d_globalData, sizeof(My_CUDA_Global)));
	CUDA_SAFE_CALL(cudaMemcpy(d_globalData, globalData, sizeof(My_CUDA_Global), cudaMemcpyHostToDevice));	
	//CUDA_SAFE_CALL(cudaMemcpyToSymbol(d_globalData, globalData, sizeof(My_CUDA_Global)));	

	CUDA_SAFE_CALL(cudaMalloc((void**)&d_silData, sizeof(My_CUDA_Sil)));
	CUDA_SAFE_CALL(cudaMemcpy(d_silData, silData, sizeof(My_CUDA_Sil), cudaMemcpyHostToDevice));
	//CUDA_SAFE_CALL(cudaMemcpyToSymbol(d_silData, silData, sizeof(My_CUDA_Sil)));

	CUDA_SAFE_CALL(cudaMalloc((void**)&d_silX, dSize*sizeof(int)));
	CUDA_SAFE_CALL(cudaMemcpy(d_silX, silData->silX, dSize*sizeof(int), cudaMemcpyHostToDevice));
	//CUDA_SAFE_CALL(cudaMemcpyToSymbol(d_silX, silData->silX, dSize*sizeof(int)));
	

	CUDA_SAFE_CALL(cudaMalloc((void**)&d_ZStep, dSize/2*sizeof(float)));
	//CUDA_SAFE_CALL(cudaMemcpy(d_ZStep, silData->ZStep, dSize/2*sizeof(float), cudaMemcpyHostToDevice));
	CUDA_SAFE_CALL(cudaMemcpyToSymbol(d_ZStep, silData->ZStep, dSize/2*sizeof(float)));

	CUDA_SAFE_CALL(cudaMalloc((void**)&d_silZ, dSize*sizeof(float)));
	CUDA_SAFE_CALL(cudaMemcpy(d_silZ, silData->silZ, dSize*sizeof(float), cudaMemcpyHostToDevice));
	//CUDA_SAFE_CALL(cudaMemcpyToSymbol(d_silZ, silData->silZ, dSize*sizeof(float)));

	printf("dSize = %d\n", dSize);
	CUDA_SAFE_CALL(cudaMalloc((void**)&d_silUVW, dSize*3*sizeof(float)));
	//CUDA_SAFE_CALL(cudaMemcpy(d_silUVW, silData->silUVW, dSize*3*sizeof(float), cudaMemcpyHostToDevice)); 
	CUDA_SAFE_CALL(cudaMemcpyToSymbol(d_silUVW, silData->silUVW, dSize*3*sizeof(float))); 
}

extern "C"
void WriteBackData(My_CUDA_Sil* silData, My_CUDA_Global* globalData)
{
    cudaMemcpy(globalData, d_globalData, sizeof(My_CUDA_Global), cudaMemcpyDeviceToHost);
    cudaMemcpy(silData, d_silData, sizeof(My_CUDA_Sil), cudaMemcpyDeviceToHost);
    int dSize=silData->silN*(silData->maxL-silData->minL);
	cudaMemcpy(silData->silX, d_silX, dSize*sizeof(int), cudaMemcpyDeviceToHost);

	cudaMemcpy(silData->ZStep, d_ZStep, dSize/2*sizeof(float), cudaMemcpyDeviceToHost);

	cudaMemcpy(silData->silZ, d_silZ, dSize*sizeof(float), cudaMemcpyDeviceToHost);

	cudaMemcpy(silData->silUVW, d_silUVW, dSize*3*sizeof(float), cudaMemcpyDeviceToHost); 
}

extern "C"
void CloseSilRender()
{
	( cudaFree(d_silX) );
	( cudaFree(d_ZStep) );
	( cudaFree(d_silZ) );
	( cudaFree(d_silUVW) );
	( cudaFree(d_globalData) );
	( cudaFree(d_silData) );
}
